use "esg_gexport_ex_to_im.dta", clear
gen tradeflow = string(excountry) + "_" + string(imcountry)
egen flow_id = group(tradeflow)
sort flow_id
gen esg_gap = ex_esg_score - im_esg_score
gen esg_gap2 = esg_gap * esg_gap
gen esg_gap3 = esg_gap2 * esg_gap
gen abs_gap = abs(esg_gap)
gen abs_gap3 = abs_gap * abs_gap * abs_gap
gen var3 = ex_esg_score * esg_gap
gen var4 = ex_esg_score * esg_gap2
gen var5 = ex_esg_score * abs_gap
drop if im_esg_score == .

duplicates drop flow_id year, force
xtset flow_id year
scatterfit g_export_rate esg_gap if esg_gap > 0, binned nq(500)
scatterfit g_export_rate esg_gap if esg_gap < 0, binned nq(500)

*Benchmark regression
reghdfe g_export_rate esg_gap if esg_gap > 0, a(excountry imcountry year) cluster(flow_id)
est store benchmark_1
reghdfe g_export_rate esg_gap if esg_gap < 0, a(excountry imcountry year) cluster(flow_id)
est store benchmark_5

rename excode code
merge m:m code year using "controls.dta", keepusing(CO2_emissions GDP_percapitaPPP urban_population exchange_rate merchandise_trade)
keep if _merge == 3
drop _merge
rename CO2_emissions CO2_x
rename GDP_percapitaPPP GDP_x
rename urban_population Urban_x
rename exchange_rate Exchange_x
rename merchandise_trade Merchan_x
rename code excode
rename imcode code
merge m:m code year using "controls.dta", keepusing(CO2_emissions GDP_percapitaPPP urban_population exchange_rate merchandise_trade)
keep if _merge == 3
drop _merge
rename CO2_emissions CO2_m
rename GDP_percapitaPPP GDP_m
rename urban_population Urban_m
rename exchange_rate Exchange_m
rename merchandise_trade Merchan_m
replace GDP_x = log(GDP_x)
replace GDP_m = log(GDP_m)
replace Exchange_m = log(Exchange_m)
replace Exchange_x = log(Exchange_x)
replace CO2_x = log(CO2_x)
replace CO2_m = log(CO2_m)
save "reg_data.dta", replace

reg g_export_rate esg_gap GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap > 0
asdoc vif, save(vif.doc) replace
reghdfe g_export_rate esg_gap GDP_x CO2_x Urban_x Exchange_x Merchan_x if esg_gap > 0, a(excountry imcountry year) vce(r)
est store benchmark_2
reghdfe g_export_rate esg_gap GDP_x CO2_x Urban_x Exchange_x Merchan_x if esg_gap < 0, a(excountry imcountry year) vce(r)
est store benchmark_6
reghdfe g_export_rate esg_gap GDP_m CO2_m Urban_m Exchange_m Merchan_m if esg_gap > 0, a(excountry imcountry year) vce(r)
est store benchmark_3
reghdfe g_export_rate esg_gap GDP_m CO2_m Urban_m Exchange_m Merchan_m if esg_gap < 0, a(excountry imcountry year) vce(r)
est store benchmark_7
reghdfe g_export_rate esg_gap GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap > 0, a(excountry imcountry year) vce(r)
est store benchmark_4
reghdfe g_export_rate esg_gap GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap < 0, a(excountry imcountry year) vce(r)
est store benchmark_8
*PPML
gen zero_flag = (g_export_rate == 0)
summarize zero_flag if esg_gap > 0
summarize zero_flag if esg_gap < 0
ppmlhdfe g_export_rate esg_gap GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap > 0, a(excountry imcountry year) vce(r)
est store ppml_1

esttab benchmark_1 benchmark_2 benchmark_3 benchmark_4 ppml_1 benchmark_5 benchmark_6 benchmark_7 benchmark_8 using benchmark.rtf, replace b(3) t starlevels(* 0.1 ** 0.05 *** 0.01)

*roubustness test
reghdfe g_export_rate esg_gap GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap > 0, a(year flow_id) vce(r)
est store est1

rename code imcode
merge m:m excode imcode year using "esg_gexport_im_from_ex.dta", keepusing(g_import_rate Gimport tradeflow)
drop if _merge == 2
drop _merge
gen NERCA = g_export_rate - g_import_rate
reghdfe NERCA esg_gap GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap > 0, a(flow_id year) vce(r)
est store est2

gen G7 = inlist(imcode, "CAN", "FRA", "DEU", "ITA", "JPN", "GBR", "USA")
rename G7 G7_m
gen G7_x = inlist(excode, "CAN", "FRA", "DEU", "ITA", "JPN", "GBR", "USA")
reghdfe g_export_rate esg_gap GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap > 0 & G7_x == 0, a(year flow_id) vce(r)
est store est3
esttab est1 est2 est3 using robustness.rtf, replace b(3) t starlevels(* 0.1 ** 0.05 *** 0.01)

*endogeneity
use "reg_data.dta", clear
merge m:1 code year using "forest_rate.dta", keepusing(Forest_rate)
drop if _merge == 2
drop _merge
rename Forest_rate Forest_rate_m
rename code imcode
rename excode code
merge m:1 code year using "forest_rate.dta", keepusing(Forest_rate)
drop if _merge == 2
drop _merge
rename Forest_rate Forest_rate_x
gen Forestgap = Forest_rate_x - Forest_rate_m
reghdfe esg_gap Forestgap if esg_gap > 0, a(excountry imcountry year)
ivreghdfe g_export_rate (esg_gap = Forestgap) GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap > 0, absorb(flow_id year) cluster(flow_id) first
est store iv
esttab iv using iv.rtf, replace b(3) t starlevels(* 0.1 ** 0.05 *** 0.01)

*mechanism
use "reg_data.dta", clear
merge m:1 code year using "oecd_eps.dta", keepusing(eps)
drop if _merge == 2
drop _merge
rename eps eps_m
gen var6 = eps_m * esg_gap
rename code imcode
rename excode code
merge m:1 code year using "oecd_eps.dta", keepusing(eps)
drop if _merge == 2
drop _merge
rename eps eps_x
gen var7 = eps_x * esg_gap
reghdfe g_export_rate esg_gap var6 eps_m GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x  Exchange_m Merchan_x Merchan_m if esg_gap > 0, a(excountry imcountry year) vce(r)
est store eps1
reghdfe g_export_rate esg_gap var7 eps_x GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap > 0, a(excountry imcountry year) vce(r)
est store eps2
esttab eps1 eps2 using eps.rtf, replace b(3) t starlevels(* 0.1 ** 0.05 *** 0.01)

merge m:1 code year using "environmental_patent.dta", keepusing(patent)
drop if _merge == 2
drop _merge
rename patent patent_x
rename code excode
rename imcode code
merge m:1 code year using "environmental_patent.dta", keepusing(patent)
drop if _merge == 2
drop _merge
rename patent patent_m
replace patent_x = log(patent_x + 1)
replace patent_m = log(patent_m + 1)
gen GI_gap = patent_x - patent_m
reghdfe GI_gap esg_gap GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap > 0, a(excountry imcountry year) vce(r)
est store GI1
reghdfe g_export_rate GI_gap esg_gap GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap > 0, a(excountry imcountry year) vce(r)
est store GI2
esttab GI1 GI2 using GI_gap.rtf, replace b(3) t starlevels(* 0.1 ** 0.05 *** 0.01)

*non-linear analysis
reghdfe g_export_rate abs_gap esg_gap2 GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap < 0, a(excountry imcountry year) vce(r)
est store nonlinear1
reghdfe g_export_rate abs_gap esg_gap2 abs_gap3 GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap < 0, a(excountry imcountry year) vce(r)
est store nonlinear2
reghdfe g_export_rate abs_gap esg_gap2 GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap < 0 & esg_gap > -0.1276, a(excountry imcountry year) vce(r)
est store nonlinear3
reghdfe g_export_rate abs_gap esg_gap2 abs_gap3 GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap < 0 & esg_gap > -0.1276, a(excountry imcountry year) vce(r)
est store nonlinear4
ppmlhdfe g_export_rate abs_gap esg_gap2 GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap < 0 & esg_gap > -0.1276, a(excountry imcountry year) vce(r)
est store ppml_2
esttab nonlinear1 nonlinear2 nonlinear3 ppml_2 nonlinear4 using nonlinear.rtf, replace b(3) t starlevels(* 0.1 ** 0.05 *** 0.01)

reghdfe g_export_rate esg_gap esg_gap2 GDP_x GDP_m CO2_x CO2_m Urban_x Urban_m Exchange_x Exchange_m Merchan_x Merchan_m if esg_gap < 0 & esg_gap > -0.1276, a(excountry imcountry year) vce(r)
utest esg_gap esg_gap2, fieller min(-0.1276) max(0) level(95)

twoway (function y = -20.69661 * (abs(x))^3 + 5.239883 * (abs(x))^2 - 0.326368 * (abs(x)), range(-.2092437 0) ) ///
, xline(-0.2092437, lpattern(dash) lcolor(gs10)) ///
xline(-0.1276, lpattern(dash) lcolor(gs10)) ///
xline(-0.0412, lpattern(dash) lcolor(gs10)) ///
xline(0, lpattern(dash) lcolor(gs10)) ///
xlabel(-0.2092437 "-0.2092" -0.1276 "-0.1276" -0.0412 "-0.0412" 0 "0", labsize(small) angle(45)) //

twoway (function y = 1.72229 * x^2 + 0.1853101 * x, range(-.1276 0) ), ///
xlabel( -0.1276 "-0.1276" -0.0537976 "-0.0538" 0 "0", labsize(small) angle(45)) //
